% rundtc1cphmon
%

% clear and close if memory, display are issues
%clear all;
%close all;

recon.Ncpus = 1; % Number of CPUs available for parallel computing, '1' for serial computing
coex=-2i*pi;
% characterize imaging experiment
pseq.sw = 198511.166253;  % sampling frequency
pseq.fovcm = 12.8;  % field of view in cm?
pseq.dcoff = 22.4+38.7i;  % DC offset in FID
pseq.offr = 35;  % off-resonance frequency in Hz
pseq.startx = 230;  % ignore this many points at start of FID
pseq.endx = 763;  % ignore this many points at end of FID

% set up image model structure
reso = 72 ;  % image model side dimension

wfi = 1.0;  % relative scaling of exponential imaginary
wfr = 1.0;  % relative scaling of exponential real
model = model_setup(reso,pseq.fovcm,wfr+1i*wfi);
% fid_file = './Dfiles/raw_data/s_20071015_01/data/pinecone_01';
% trajectory_file = 'KTR100735';
fid_file = '/home/twieg/matfeb/PARSE/Dfiles/raw_data/kiwi_pinecone_functional_saccade_s_20070112_02/data/pinecone_01';
trajectory_file = 'KTR070125';
% calibration_file = '';  % calibration data in trajectory_file

% parameters for reconstruction algorithm
recon.alphaM = 0;  % regularization parameter for magnitude
recon.alphaf = 0;  % regularization parameter for decay/frequency map
% complex alphaf applies real part to decay and imaginary to freq
recon.wfi = model.wfi;  % preconditioning weight for frequency (?)
recon.wfr = model.wfr;  % preconditioning weight for decay (?)
recon.NLIST = 120*[1;2;3;5;7;10;17;26;35;45;55;75;100];  % data lengths in progressive-length CG (PLCG)
%recon.FLIST = [.1;.08;.07;.06;.06;.05;.05;.05;.04;.03;.02;.02;.01]; 
recon.FLIST = [.1;.08;.07;.05;.04;.04;.04;.04;.04;.03;.02;.02;.01];
recon.NLENGTHS = size(recon.NLIST,1);  % # of outer loops in PLCG
%recon.FLIST = .02*ones(recon.NLENGTHS,1);
recon.NIT = 120;  % # of iterations
recon.GSS = 0;  % use golden section for line search
recon.del = 4e-09;  % step size at which to measure gradient in Cost function CC
%recon.bounddecay = 5000/pseq.sw;  % threshold to bound decay parameter
recon.magdown = 0.05;  % reduce mag by this factor if decay exceeded
recon.plot = 0;  % show plots on the fly every recon.plot iterations (0=off)
recon.vis = 2;  % show images on the fly every recon.vis iterations (0=off)
recon.val = 0;  % show criterion values on the fly
recon.time = 0;  % show times on the fly
recon.mat = 'matrix';  % matrix, separable, or filter; save computation or memory

% load data and/or calculate phantom FID and trajectory
[pseq.FR,pseq.FI] = LOAD_FID(fid_file);
eval(['load ' trajectory_file]);  % trajectory data and eddy current phases
% phip and phir computed from calibration data using krdpenu19
pseq.krd = real(kss);  % readout k-coordinate
pseq.kpe = imag(kss);  % phase-encode k-coordinate
pseq.phir = phi;  % measured eddy current phase -- r
pseq.phip = 0;  % measured eddy current phase -- p
kr=abs(kss);

% set up useful parameters and preprocess data for reconstruction
ffac=1.1*model.N/(4*pi);
frno=3;
fid = fid_setupss(pseq,recon.Ncpus,ffac,frno);
NF=size(fid.sig,2);
compphan = 1;
%=======================================================
if compphan == 1
%this section replaces experimental fid with computed "clock" phantom
Rout=5;Rin=1;radi=4;
%xoff=4;yoff=1;
froff1=0;
alphalist=(0:11)'*pi/6 + pi/12;
frofflist=[zeros(6,1);100*ones(6,1)];
amplist=[0.5;0.2;0.1;0.05;0.02;0.01;0.01;0.02;0.05;0.1;0.2;0.5];
FIDctr=(Rout^2.*besselj(1,pi*Rout.*kr)./(pi*Rout.*kr))...
   .*exp((-30+2i*pi*froff1)*(1:NF)'*fid.delt);
FRIS=FIDctr;
for nclock=1:12
    alp=alphalist(nclock);xoff=radi*cos(alp);yoff=radi*sin(alp);
froff2=frofflist(nclock);amp=amplist(nclock);
FRIS=FRIS+amp*exp(2i*pi*(pseq.krd*xoff+pseq.kpe*yoff+froff2*(1:NF)'*fid.delt))...
    .*(Rin^2.*besselj(1,pi*Rin.*kr)./(pi*Rin.*kr)).*exp((-30)*(1:NF)'*fid.delt);
end
pseq.FR=real(FRIS);pseq.FI=imag(FRIS);
fid = fid_setupscph(pseq,recon.Ncpus,ffac);
noisefac=.000001;
fid.sig = fid.sig + (noisefac*max(abs(fid.sig)))*(randn(1,NF)+1i*randn(1,NF));
pseq.FR=real(fid.sig);pseq.FI=imag(fid.sig);
end
%===========================================
%figure(7);plot(real(fid.sig));
%Nstages=22;%XXX -- latest monkey trial!
Nstages=20;
M0Cc=zeros(reso,reso,Nstages);R2ec=M0Cc;frmapc=M0Cc;fidrescale=zeros(1,Nstages);
%M0Cd=M0Cc;R2ed=M0Cc;frmapd=M0Cc;
NLL=length(model.ix);
kbasis = make_kbasis(fid,model,'matrix',recon.mat,recon.Ncpus);
model.pvec = zeros(1,2*NLL);
modelsav = model;
NM = model.N;
%fraclist=0.6*ones(1,Nstages);
%fraclist=0.4+0.5*(0:Nstages-1)/(Nstages-1);%XXX -- latest monkey trial!
fraclist=0.5+0.4*(0:Nstages-1)/(Nstages-1);%XXX -- latest monkey trial!
for nrun=1:Nstages
    modelsav.pvec = zeros(1,2*NLL);
    %modelsav.pvec = [0.0001*ones(1,NLL) zeros(1,NLL)];
    model = cgparse5(fid,recon,modelsav,kbasis);
    %model.pvec(NM+1:end)=-min(0,real(model.pvec(NM+1:end))) + 1i*imag(model.pvec(NM+1:end));
    vec_ind = model.ix+model.reso*(model.iy-1);
    model.amp = zeros(model.reso);
    model.exp = zeros(model.reso);
    model.amp(vec_ind) = model.pvec(1:NM);
    model.exp(vec_ind) = model.pvec(NM+1:end);
    model2=model;
    frac=fraclist(nrun);
    MTH1 = (1-exp(wfr*real(model.pvec((NM+1):end)).*NF)).*abs(model.pvec(1:NM))./abs(real(model.pvec((NM+1):end))+eps);
    MTH=MTH1./max(MTH1);%MTH represents area under the amplitude of the local signal curve
    %here model magnitudes are scaled by (1-frac) to represent portion
    %subtracted from fid below
    model2.pvec(1:NM)=(1-frac)*(MTH>frac).*model.pvec(1:NM);
    model2.amp(vec_ind) = model2.pvec(1:NM);
    model2.exp(vec_ind) = model2.pvec(NM+1:end);
    %subtract contribution of "top part" of object from the fid:
    fidnext = fid.sig - fidsyn(model2,kbasis,length(fid.sig));
    %fid.sig = fidnext;
     fidrescale(nrun)=1*ffac./sum(abs(fidnext));
     fid.sig = fidrescale(nrun)*fidnext;
    
    pseq.FR=real(fid.sig);pseq.FI=imag(fid.sig);
     model2.amp = zeros(model.reso);
     model2.exp = zeros(model.reso);
%      model.amp = zeros(model.reso);
%      model.exp = zeros(model.reso);
    model2.amp(vec_ind) = model2.pvec(1:NM);
    model2.exp(vec_ind) = model2.pvec(NM+1:end);
%     model.amp(vec_ind) = model.pvec(1:NM);
%     model.exp(vec_ind) = model.pvec(NM+1:end);
    %----------------------
    R2ec(:,:,nrun) = real(model2.exp)*(-model.wfr/fid.delt);
    frmapc(:,:,nrun) = imag(model2.exp)*(model.wfi/(2*pi*fid.delt));
    M0Cc(:,:,nrun) = model2.amp;
%     R2ed(:,:,nrun) = real(model.exp)*(-model.wfr/fid.delt);
%     frmapd(:,:,nrun) = imag(model.exp)*(model.wfi/(2*pi*fid.delt));
%     M0Cd(:,:,nrun) = model.amp;
recon.NLIST = 120*[1;2;3;5;7;9;12;15;20;28;36;44;55;75;100];  % data lengths in progressive-length CG (PLCG)
recon.FLIST = [.1;.09;.07;.05;.04;.04;.03;.03;.03;.03;.03;.03;.03;.02;.01];
recon.NLENGTHS = size(recon.NLIST,1);
%modelsav.pvec=model2.pvec;
    %model.pvec = 0.01*model2.pvec;%added 12/21 late
end
%figure(8);for m=1:15;subplot(3,5,m);imagesc(abs(M0Cc(:,:,m)));axis off;end

%process for Mest, R2est, frest:
Mest=zeros(model.reso);frest=Mest;R2est=frest;netrescale=1;

for m=1:Nstages
    Mest=Mest+(Mest==0).*M0Cc(:,:,m)/netrescale;
    frest=frest+(frest==0).*(Mest~=0).*frmapc(:,:,m);
    R2est=R2est+(R2est==0).*(Mest~=0).*R2ec(:,:,m);
    netrescale=netrescale*fidrescale(m);
end
Mest=Mest./max(abs(Mest(:)));
figure(10);subplot(1,3,1);imagesc(abs(Mest));axis image;axis off;
subplot(1,3,2);imagesc(R2est);axis image;axis off;
subplot(1,3,3);imagesc(frest);axis image;axis off;